In [487]:
import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
import seaborn as sns #for plotting
from sklearn.ensemble import RandomForestClassifier #for the model
from sklearn.tree import DecisionTreeClassifier
from sklearn.tree import export_graphviz #plot tree
from sklearn.metrics import roc_curve, auc #for model evaluation
from sklearn.metrics import classification_report #for model evaluation
from sklearn.metrics import confusion_matrix #for model evaluation
from sklearn.model_selection import train_test_split #for data splitting
import eli5 #for purmutation importance
from eli5.sklearn import PermutationImportance
from sklearn.linear_model import LogisticRegression
from sklearn.ensemble import RandomForestClassifier
from sklearn.neural_network import MLPClassifier
from sklearn.preprocessing import LabelEncoder
from sklearn.model_selection import train_test_split
from sklearn.metrics import f1_score,accuracy_score,roc_auc_score
from xgboost import XGBClassifier, plot_importance
from pdpbox import pdp, get_dataset, info_plots
# model explainers
import lime
from lime.lime_tabular import LimeTabularExplainer
import shap
from shap import TreeExplainer,KernelExplainer,LinearExplainer
shap.initjs()
import eli5
import shap
In [488]:
pip install PDPbox
In [489]:
np.random.seed(123) #ensure reproducibility
pd.options.mode.chained_assignment = None #hide any pandas warnings
In [490]:
from pdpbox import pdp, info_plots #for partial plots
In [491]:
data = pd.read_csv("/Users/odianosenakhibi/Downloads/Thesis/data/heart.csv")
In [492]:
data.head(10)
Out[492]:
In [493]:
data.describe()
Out[493]:
In [494]:
data.info()
In [495]:
data.dtypes
Out[495]:
In [536]:
plt.figure(figsize=(10,8))
sns.heatmap(data.isnull().T, cbar=False)
Out[536]:
In [496]:
data = pd.get_dummies(data, drop_first=True)
In [497]:
data.head()
Out[497]:
In [498]:
sns.countplot(x='target',data=data)
Out[498]:
In [499]:
sns.pairplot(data)
plt.show()
In [500]:
# Show the results of a linear regression within each dataset
sns.lmplot(x="trestbps", y="chol",data=data,hue="cp")
plt.show()
In [501]:
sns.violinplot(data.age, palette="Set3", bw = .2, cut = 1, linewidth = 1)
plt.xticks(rotation = 90)
plt.title("Age Rates")
plt.show()
In [502]:
plt.figure(figsize=(15,7))
sns.violinplot(x=data.age,y=data.target)
plt.xticks(rotation=90)
plt.legend()
plt.title("Age & Target System")
plt.show()
In [503]:
#split the data
X_train, X_test, y_train, y_test = train_test_split(data.drop('target', 1), data['target'], test_size = .2, random_state=10)
In [504]:
#Let's see how the correlation values between them
data.corr()
Out[504]:
In [505]:
model = RandomForestClassifier(max_depth = 5)
model.fit(X_train, y_train)
Out[505]:
In [506]:
print('X_train',X_train.shape)
print('X_test',X_test.shape)
print('y_train',y_train.shape)
print('y_test',y_test.shape)
In [507]:
#Normalization
X_train=(X_train-np.min(X_train))/(np.max(X_train)-np.min(X_train)).values
X_test=(X_test-np.min(X_test))/(np.max(X_test)-np.min(X_test)).values
In [508]:
from sklearn.naive_bayes import GaussianNB
from sklearn.neighbors import KNeighborsClassifier
from sklearn.tree import DecisionTreeClassifier
from sklearn.ensemble import RandomForestClassifier
from sklearn.ensemble import GradientBoostingClassifier
from sklearn.svm import SVC
from sklearn.metrics import roc_curve, auc
from sklearn.model_selection import GridSearchCV,train_test_split,cross_val_score
import os
from sklearn.decomposition import PCA
pca=PCA().fit(X_train)
print(pca.explained_variance_ratio_)
print()
print(X_train.columns.values.tolist())
print(pca.components_)
In [509]:
#To reduce the data dimensions and noise in the input data
pca = PCA(n_components=8)
pca.fit(X_train)
reduced_data_train = pca.transform(X_train)
#inverse_data = pca.inverse_transform(reduced_data)
plt.scatter(reduced_data_train[:, 0], reduced_data_train[:, 1],
label='reduced')
plt.xlabel('First Principal Component')
plt.ylabel('Second Principal Component')
plt.show()
In [510]:
pca = PCA(n_components=8)
pca.fit(X_test)
reduced_data_test = pca.transform(X_test)
#inverse_data = pca.inverse_transform(reduced_data)
plt.scatter(reduced_data_test[:, 0], reduced_data_test[:, 1],
label='reduced')
plt.xlabel('First Principal Component')
plt.ylabel('Second Principal Component')
plt.show()
In [511]:
estimator = model.estimators_[1]
feature_names = [i for i in X_train.columns]
y_train_str = y_train.astype('str')
y_train_str[y_train_str == '0'] = 'no disease'
y_train_str[y_train_str == '1'] = 'disease'
y_train_str = y_train_str.values
In [512]:
export_graphviz(estimator, out_file='tree.dot',
feature_names = feature_names,
class_names = y_train_str,
rounded = True, proportion = True,
label='root',
precision = 2, filled = True)
from subprocess import call
call(['dot', '-Tpng', 'tree.dot', '-o', 'tree.png', '-Gdpi=600'])
from IPython.display import Image
Image(filename = 'tree.png')
#to capture the interaction between the features of the data
Out[512]:
In [513]:
y_predict = model.predict(X_test)
y_pred_quant = model.predict_proba(X_test)[:, 1]
y_pred_bin = model.predict(X_test)
In [514]:
confusion_matrix = confusion_matrix(y_test, y_pred_bin)
confusion_matrix
Out[514]:
In [515]:
total = sum(sum(confusion_matrix))
sensitivity = confusion_matrix[0, 0] / (confusion_matrix[0, 0] + confusion_matrix[1, 0])
print('Sensitivity : ', sensitivity)
specificity = confusion_matrix[1, 1] / (confusion_matrix[1, 1] + confusion_matrix[0, 1])
print('Specificity : ', specificity)
In [516]:
fpr, tpr, thresholds = roc_curve(y_test, y_pred_quant)
fig, ax = plt.subplots()
ax.plot(fpr, tpr)
ax.plot([0, 1], [0, 1], transform=ax.transAxes, ls="--", c=".3")
plt.xlim([0.0, 1.0])
plt.ylim([0.0, 1.0])
plt.rcParams['font.size'] = 12
plt.title('ROC curve for Heart Disease classifier')
plt.xlabel('False Positive Rate (1 - Specificity)')
plt.ylabel('True Positive Rate (Sensitivity)')
plt.grid(True)
In [517]:
auc(fpr, tpr)
Out[517]:
In [518]:
import seaborn as sns
In [519]:
from xgboost import XGBClassifier, plot_importance
In [520]:
import warnings
warnings.filterwarnings('ignore')
plt.style.use('fivethirtyeight')
%matplotlib inline
In [521]:
# for visualizing correlations
f, ax = plt.subplots(figsize=(10, 6))
corr = data.corr()
hm = sns.heatmap(round(corr,2), annot=True, ax=ax, cmap="Reds",fmt='.2f',
linewidths=.05)
f.subplots_adjust(top=0.93)
t= f.suptitle('Heart Disease Attributes Correlation Heatmap', fontsize=14)
In [522]:
target = 'target'
features_list = list(data.columns)
features_list.remove(target)
In [523]:
sns.set()
sns.relplot(data = data, x = 'cp', y = target, kind = 'line', height = 5, aspect = 2, color = 'red')
Out[523]:
In [524]:
%%time
# ML in two lines ;)
xgb = XGBClassifier(objective='binary:logistic', random_state=33, n_jobs=-1)
xgb.fit(X_train, y_train)
In [525]:
xgb = XGBClassifier(objective='binary:logistic', random_state=33, n_jobs=-1)
xgb.fit(X_train, y_train)
Out[525]:
In [526]:
xgb_predictions = xgb.predict(X_test)
In [527]:
# We design a simple classification evaluative function
def evaluation_scores(test, prediction, target_names=None):
print('Accuracy:', np.round(metrics.accuracy_score(test, prediction), 10))
print('-'*60)
print('classification report:\n\n', metrics.classification_report(y_true=test, y_pred=prediction, target_names=target_names))
classes = [0, 1]
total_classes = len(classes)
level_labels = [total_classes*[0], list(range(total_classes))]
cm = metrics.confusion_matrix(y_true=test, y_pred=prediction, labels=classes)
cm_frame = pd.DataFrame(data=cm, columns=pd.MultiIndex(levels=[['Predicted:'], classes], labels=level_labels),
index=pd.MultiIndex(levels=[['Actual:'], classes], labels=level_labels))
print('-'*60)
print('Confusion matrix:\n')
print(cm_frame)
In [528]:
from sklearn import metrics
from sklearn.metrics import accuracy_score, roc_curve, auc
In [529]:
evaluation_scores(y_test, xgb_predictions, target_names=['no disease', 'disease'])
#Showing the main classification metrics, (precision, recall and f1 score)
In [303]:
# calculate the FPR and TPR for all thresholds of the classification
probs = xgb.predict_proba(X_test)
preds = probs[:,1]
fpr, tpr, threshold = metrics.roc_curve(y_test, preds)
roc_auc = metrics.auc(fpr, tpr)
#AUC-ROC measures the performance of the classification problem at various thresholds with the ROC curve telling us the ability of the model to distinguish between classes
plt.title('Receiver Operating Characteristic')
plt.plot(fpr, tpr, 'red', label = 'ROC AUC score = %0.2f' % roc_auc)
plt.legend(loc = 'lower right')
plt.plot([0, 1], [0, 1],'b--')
plt.xlim([0, 1])
plt.ylim([0, 1])
plt.ylabel('True Positive Rate')
plt.xlabel('False Positive Rate')
plt.show()
In [535]:
auc(fpr,tpr)
Out[535]:
In [ ]:
In [304]:
# ploting XGBoost default feature importances
fig = plt.figure(figsize = (18, 10))
title = fig.suptitle("Native Feature Importances from XGBoost", fontsize=14)
ax1 = fig.add_subplot(2, 2, 1)
plot_importance(xgb, importance_type='weight', ax=ax1, color='red')
ax1.set_title("Feature Importance with Feature Weight");
ax2 = fig.add_subplot(2, 2, 2)
plot_importance(xgb, importance_type='cover', ax=ax2, color='red')
ax2.set_title("Feature Importance with Sample Coverage");
ax3 = fig.add_subplot(2, 2, 3)
plot_importance(xgb, importance_type='gain', ax=ax3, color='red')
ax3.set_title("Feature Importance with Split Mean Gain");
In [305]:
import eli5
from eli5.sklearn import PermutationImportance
In [306]:
eli5.show_weights(xgb.get_booster())
#this gives us an idea of how the model performs globally
Out[306]:
In [58]:
%%time
# we need to retrain a new model with arrays
# as eli5 has a bug with Dataframes and XGBoost
# cf. https://github.com/TeamHG-Memex/eli5/pull/261
xgb_array = XGBClassifier(objective='binary:logistic', random_state=33, n_jobs=-1)
xgb_array.fit(X_train, y_train)
In [59]:
xgb_array = XGBClassifier(objective='binary:logistic', random_state=33, n_jobs=-1)
xgb_array.fit(X_train, y_train)
Out[59]:
In [219]:
feat_permut = PermutationImportance(xgb_array,
random_state=1).fit(X_train, y_train)
eli5.show_weights(feat_permut, feature_names = features_list)
Out[219]:
In [62]:
from pdpbox import pdp, get_dataset, info_plots
In [63]:
def plot_pdp(model, df, feature, cluster_flag=False, nb_clusters=None, lines_flag=False):
# Create the data that we will plot
pdp_goals = pdp.pdp_isolate(model=model, dataset=data, model_features=df.columns.tolist(), feature=feature)
# plot it
pdp.pdp_plot(pdp_goals, feature, cluster=cluster_flag, n_cluster_centers=nb_clusters, plot_lines=lines_flag)
plt.show()
In [64]:
# plot the PD univariate plot
plot_pdp(xgb, X_train, 'cp')
In [65]:
# for ICE plot we must specify the numbers of similarity clusters we want
# here 24
plot_pdp(xgb, X_train, 'cp', cluster_flag=True, nb_clusters=24, lines_flag=True)
In [68]:
features_to_plot = ['thalach', 'age']
inter1 = pdp.pdp_interact(model=xgb, dataset=X_train, model_features=features_list, features=features_to_plot)
pdp.pdp_interact_plot(pdp_interact_out=inter1, feature_names=features_to_plot, plot_type='grid')
# we use plot_type='grid' as the default and better option 'contour' has a bug which is being corrected
# cf. https://github.com/SauceCat/PDPbox/issues/40
plt.show()
In [83]:
cat_cols = list(data.select_dtypes('object').columns)
class_dict = {}
for col in cat_cols:
data = pd.concat([data.drop(col, axis=1), pd.get_dummies(data[col])], axis=1)
data.head()
Out[83]:
In [84]:
ML_models = {}
model_index = ['LR','RF','NN']
model_sklearn = [LogisticRegression(solver='liblinear',random_state=0),
RandomForestClassifier(n_estimators=100,random_state=0),
MLPClassifier([100]*5,early_stopping=True,learning_rate='adaptive',random_state=0)]
model_summary = []
for name,model in zip(model_index,model_sklearn):
ML_models[name] = model.fit(X_train,y_train)
preds = model.predict(X_test)
model_summary.append([name,f1_score(y_test,preds,average='weighted'),accuracy_score(y_test,preds),
roc_auc_score(y_test,model.predict_proba(X_test)[:,1])])
ML_models
Out[84]:
In [85]:
model_summary = pd.DataFrame(model_summary,columns=['Name','F1_score','Accuracy','AUC_ROC'])
model_summary = model_summary.reset_index()
display(model_summary)
In [86]:
g=sns.regplot(data=model_summary, x="index", y="AUC_ROC", fit_reg=False,
marker="o", color="black", scatter_kws={'s':500})
for i in range(0,model_summary.shape[0]):
g.text(model_summary.loc[i,'index'], model_summary.loc[i,'AUC_ROC']+0.02, model_summary.loc[i,'Name'],
horizontalalignment='center',verticalalignment='top', size='large', color='black')
In [482]:
#initialization of an explainer from LIME
explainer = LimeTabularExplainer(X_train.values,
mode='classification',
feature_names=X_train.columns,
class_names=['disease','no disease'])
In [483]:
exp = explainer.explain_instance(X_test.head(1).values[0],
ML_models['LR'].predict_proba,
num_features=X_train.shape[1])
exp.show_in_notebook(show_table=True, show_all=True)
In [484]:
exp = explainer.explain_instance(X_test.head(1).values[0],
ML_models['RF'].predict_proba,
num_features=X_train.shape[1])
exp.show_in_notebook(show_table=True, show_all=False)
In [485]:
exp = explainer.explain_instance(X_test.head(1).values[0],
ML_models['NN'].predict_proba,
num_features=X_train.shape[1])
exp.show_in_notebook(show_table=True, show_all=False)
In [91]:
#ELI5
eli5.show_weights(ML_models['LR'], feature_names = list(X_test.columns))
Out[91]:
In [92]:
eli5.show_prediction(ML_models['LR'], X_test.head(1).values[0],feature_names=list(X_test.columns))
Out[92]:
In [93]:
exp = PermutationImportance(ML_models['LR'],
random_state = 0).fit(X_test, y_test)
eli5.show_weights(exp,feature_names=list(X_test.columns))
Out[93]:
In [94]:
eli5.show_weights(ML_models['RF'],feature_names=list(X_test.columns))
Out[94]:
In [95]:
eli5.show_prediction(ML_models['RF'], X_test.head(1).values[0],feature_names=list(X_test.columns))
Out[95]:
In [96]:
exp = PermutationImportance(ML_models['RF'],
random_state = 0).fit(X_test, y_test)
eli5.show_weights(exp,feature_names=list(X_test.columns))
Out[96]:
In [97]:
eli5.show_weights(ML_models['NN'])
Out[97]:
In [531]:
explainer = shap.TreeExplainer(xgb)
shap_values = explainer.shap_values(X_test)
In [532]:
X_shap = pd.DataFrame(shap_values)
X_shap.tail()
Out[532]:
In [533]:
print('Expected Value: ', explainer.expected_value)
In [534]:
shap.summary_plot(shap_values, X_test, plot_type="bar", color= 'red')
In [98]:
#SHAP
explainer = LinearExplainer(ML_models['LR'], X_train, feature_dependence="independent")
shap_values = explainer.shap_values(X_test.head(1).values)
shap.force_plot(explainer.expected_value,
shap_values,
X_test.head(1).values,
feature_names=X_test.columns)
Out[98]:
In [99]:
shap_values = explainer.shap_values(X_test.head(250).values)
shap.force_plot(explainer.expected_value,
shap_values,
X_test.head(250).values,
feature_names=X_test.columns)
Out[99]:
In [100]:
shap_values = explainer.shap_values(X_test.values)
spplot = shap.summary_plot(shap_values, X_test.values, feature_names=X_test.columns)
In [101]:
top4_cols = ['ca','cp','thalach','sex']
for col in top4_cols:
shap.dependence_plot(col, shap_values, X_test)
In [102]:
explainer = TreeExplainer(ML_models['RF'])
shap_values = explainer.shap_values(X_test.head(1).values)
shap.force_plot(explainer.expected_value[1],
shap_values[1],
X_test.head(1).values,
feature_names=X_test.columns)
Out[102]:
In [104]:
X_train_kmeans = shap.kmeans(X_train, 10)
explainer = KernelExplainer(ML_models['NN'].predict_proba,X_train_kmeans)
shap_values = explainer.shap_values(X_test.head(1).values)
shap.force_plot(explainer.expected_value[1],
shap_values[1],
X_test.head(1).values,
feature_names=X_test.columns)
Out[104]:
In [ ]: